Topological phase diagrams and Majorana zero modes of the Kitaev ladder and tube
Wang Yiming1, Li Zhidan1, Han Qiang1, 2, †
Department of Physics, Renmin University of China, Beijing 100872, China
Beijing Key Laboratory of Opto-electronic Functional Materials and Micro-nano Devices, Renmin University of China, Beijing 100872, China

 

† Corresponding author. E-mail: hanqiang@ruc.edu.cn

Abstract

In this paper, we study two quasi-one-dimensional (1D) Kitaev models with ladder-like and tube-like spatial structures, respectively. Our results provide the phase diagrams and explicit expressions of the Majorana zero modes. The topological phase diagrams are obtained by decomposing the topological invariants and the topological conditions for topologically nontrivial phases are given precisely. For systems which belongs to topological class BDI, we obtain the regions in the phase diagrams where the topological numbers show even-odd effect. For the Kitaev tube model a phase factor induced by the magnetic flux in the axial direction of the tube is introduced to alter the classification of the tube Hamiltonian from class BDI to D. The Kitaev tube of class D is characterized by the index when the number of chains is odd while 0, 1, 2 when the number of chains is even. The phase diagrams show periodic behaviors with respect to the magnetic flux. The bulk-boundary correspondence is demonstrated by the observations that the topological conditions for the bulk topological invariant to take nontrivial values are precisely those for the existence of the Majorana zero modes.

1. Introduction

In 2001, Kitaev proposed a one-dimensional (1D) spinless p-wave superconductor,[1] which is the first and simplest prototype to realize the long-sought Majorana fermions[2] in solid state systems. This so-called Kitaev chain was demonstrated theoretically to host an unpaired Majorana zero mode (MZM) at each of its ends.[1,310] These MZMs are exotic quasiparticle excitations which are their own antiparticles, each of which can be considered as half of a fermion. Topological superconductors supporting MZMs are candidate platforms to realize fault tolerant topological quantum computation and storage due to their non-Abelian statistics and immunity to local perturbations.[1116] Intense research interest has been stimulated in the community of condensed matter physics since the seminal paper by Kitaev. Various 1D topological superconducting systems have been proposed theoretically[1720] to realize effective p-wave pairing in real materials with spinful fermions and several experiments[2124] reported signatures of observing the MZMs in the InSb or InAs nanowire systems based on the theoretical suggestions that the semiconductor nanowires with strong Rashba spin–orbit coupling in proximity to an s-wave superconductor can support MZMs in presence of a Zeeman field.

On the other hand there are also theoretical efforts[2531] on generalizing the Kitaev chain model to quasi-1D spinless fermion systems to provide practical route to realize MZMs in thin-film microstructures[25,32] or to understand the signatures of MZMs observed in nanowire systems[30,33] which have finite width instead of a strict 1D atomic chain. A quasi-1D -wave superconductor with strip geometry was first investigated[25] and robust Majorana ends states were found beyond the strict 1D limit. Later similar ideas based on the quasi-1D ladder model were also discussed by several authors.[2831] The so-called Kitaev ladder is actually multiple Kitaev chains placed in parallel and coupled via inter-chain hopping or pairing with each other. People found that the Kitaev ladder system undergoes multiple topological phase transitions and the number of MZMs varies as a function of system parameters.

In this paper, we give some analytic results of the quasi-1D ladder-like and pipe-like Kitaev models. These so-called Kitaev ladder and Kitaev tube are generalization of the Kitaev chain model and can be used to simulate the thin-film microstructures[25,32] and nanowire systems,[2124] which are commonly employed to implement topological superconductors experimentally. The topological phase diagrams and the MZMs are investigated analytically for these two types of quasi-1D topological superconductors. For both the Kitaev ladder belonging to class BDI and the Kitaev tube, we decompose the topological invariants of the system into sums of independent topological numbers by transformations in momentum space. The dependence of the phase diagrams on the parity of the number of chains is also examined for all the cases studied in this work. We confirm the bulk-boundary correspondence by combining the results of MZMs in real space with the topological properties in momentum space.

This paper is organized as follows. In section 2, we discuss the Kitaev ladder model. The topological invariant and phase diagrams are studied. The MZMs are investigated. In section 3, the Kitaev tube model is introduced. Both the class BDI and class D cases are discussed. Finally, We give the conclusion in section 4.

2. Multi-leg Kitaev ladder

In this section we consider the multi-leg Kitaev ladder model that is a quasi-1D model of topological superconductors. The so-called Kitaev ladder is composed of identical Kitaev chains placed in parallel and coupled with each other as shown Fig. 1. The ladder Hamiltonian is given by

where M is the number of legs of the ladder and N denotes the length of each chain. and represent the annihilation and creation operators of spinless fermion on the n-th site of the m-th leg, respectively. Here and are the intra(inter)-chain hopping integral and pairing potential, respectively. μ is the chemical potential. We consider the ladder Hamiltonian belonging to topological class BDI[34] which has time-reversal, particle-hole and chiral symmetries.[30] tx, ty, , and μ are all real variables. Without loss of generality, we set . Furthermore can also be set as positive under suitable gauge.

Fig. 1. Illustration of quasi-1D Kitaev ladder model, which is composed of M identical Kitaev chains of length N coupled with each other via nearest-neighbor inter-chain hopping ty and pairing . tx and denote the intra-chain hopping integral and pairing potential.
2.1. Topological invariant and phase diagrams

To show the topological properties of this system, we take a Fourier transformation as shown in Eq. (2) along the x direction

and the Hamiltonian (1) can be expressed as with and
where and are M × M tridiagonal matrix whose elements are
The Bogoliubov–de-Gennes (BdG) Hamiltonian Hk has chiral symmetry with the Pauli matrix and therefore is converted into a block-off-diagonal form by the following unitary transformation, and we have
Here the unitary matrix U is
with IM the M × M identity matrix. Qk in Eq. (7) is given by . The topological invariant of this quasi-1D BDI system[30,35] can be expressed as
where denotes the winding number calculated according to
with the m-th eigenvalue of the tridiagonal matrix Qk, which is written as
with . Equation (9) indicates that the topological invariant of the whole Kitaev ladder system is simply the sum of M individual topological numbers. Each as shown in Eq. (10) denotes the winding number of a closed curve formed by on the complex plane around the origin. And takes 0 or 1 according to Eqs. (10) and (11). The sum of the first three terms of the right hand side of Eq. (11) is the result of a strict 1D Kitaev chain, while the last term shows the renormalization effect resulting from the coupling among chains. In the following we will discuss two cases according to the sign of .

The first case is where the real part of Eq. (11) is renormalized and consequently takes nontrivial value 1 under the following topological condition

which is identical with the result of Ref. [30]. Hence there are totally M+1 topologically distinct phases depending on how many topological conditions (12) can be fulfilled when m ranges from 1 to M. According to Ref. [30], the topological index is Z regardless of the parity of the number of legs.

The second case is , for which the topological condition of the class-BDI system is not given by the authors of Ref. [30], although they discussed the topological condition of the class-D system. In this case, the imaginary part of Eq. (11) is renormalized and accordingly we find that the topological condition becomes

which shows that topological conditions for m and are actually identical and therefore . Equation (9) can then be further reduced to
One can see that there is no contribution from the first term on the right hand side of Eq. (14) to ν if M is even, and the topological invariant ν always changes by 2 due to the second term.

The topological phase diagrams of 4- and 5-leg Kitaev ladder models belonging to class BDI are shown in Fig. 2. The system is characterized by the Z index in the region as illustrated by Fig. 2, which is compatible with the results of Ref. [30]. However, in the region where there is an even-odd effect with respect to the number of legs of the ladder as indicated by Eq. (14). To be specific, for even M the topological number is and there are topologically distinct phases. Whereas for odd M the topological number is and there are topologically distinct phases. Accordingly, this even-odd effect leads to rich phase diagrams as illustrated in Fig. 2 for the Kitaev ladder belonging to class BDI in comparison with the simple phase diagrams of that belonging to class D in the region .[30]

Fig. 2. (color online) Phase diagrams of the 4-leg and 5-leg Kitaev ladder models. tx is chosen as energy unit. and . The integer numbers in the figure represent the Z index. In the regions where , the topological indices are 0,2,4 for M = 4 (upper panel) and 0,1,3,5 for M = 5 (lower panel).
2.2. Majorana zero modes

Here we study the presence of MZMs at the boundary of the Kitaev ladder. In the following basis of Majorana fermions,

which satisfy the algebra , , , , and , the tight-binding Hamiltonian (1) is then rewritten as
For a very long (semi-infinite) ladder with , the MZMs of the system can be written in the following form
where is the wave function of the Majorana fermion to be determined. As a Majorana fermion and thus is real. The wave function satisfies the boundary conditions
in the x direction and
in the y direction. Note that the second condition of Eq. (18) is due to the localized nature of MZMs around the left boundary of the ladder. From since γ is a zero-energy mode, we obtain the following second-order difference equations,
for , which forms one M-dimensional second-order difference equation system. We take a transformation along the y direction as follows:
and then obtain M independent equations written as,
where
with . The above results show manifestly that the Kitaev ladder can be effectively decoupled into M independent Kitaev chains each of which has a renormalized chemical potential . The MZMs in the Kitaev ladder can be found from Eq. (22) for each q separately.

For the case of , we find that the solution of Eq. (22) exists iff q satisfies the topological condition (12) (see Appendix A for detail) and correspondingly there exists an MZM according to Eq. (17) after substituting into Eq. (21)

with obviously obeying the boundary conditions (18) and (19).

For , the solution of Eq. (22) exists iff the topological condition Eq. (13) is satisfied for q. Furthermore is generally complex because may have an imaginary part. Considering that equation (13) is met simultaneously for q and and that and , we have two complex fermions and with relation according to Eqs. (21) and (24). Consequently, two Majorana fermions and are obtained pairwise. Note that there exists an exception with for odd M. In this case and one MZM is built according to Eq. (24).

From the above results of MZMs, it is demonstrated that the topological conditions for the bulk topological invariant to take nontrivial values are precisely those for the existence of boundary MZMs. Thus the topologically distinct phases are characterized by either the winding numbers obtained in momentum space or equivalently by the number of MZMs in real space. These observations obviously manifest the bulk-boundary correspondence in the Kitaev ladder system.

3. Kitaev tube

In this section we consider the Kitaev tube model as shown in Fig. 3. Similar to the ladder model, the Kitaev tube is also composed of M identical Kitaev chains that are coupled with each other. The significant difference between the tube and ladder is the boundary condition along the y direction, namely we have closed boundary condition for the former and open for the latter. The merit of the tube structure is that one can change the classification of the model by applying a magnetic field along its axial direction. The tube Hamiltonian is written as

where the meanings of the operators and parameters are the same as those of the ladder model Eq. (1) except that there are two additional parameters θ and ϕ. θ denotes the relative phase between the intra- and inter-chain hopping integrals, which results from the magnetic field applied in the axial direction of the tube as shown in Fig. 3. ϕ denotes the relative phase between the intra- and inter-chain pairing potentials. If ϕ=0 we have a -wave pairing model and if we have a -wave pairing model.

Fig. 3. (color online) Illustration of the Kitaev tube which can be viewed as the Kitaev ladder (see Fig. 1) rolled up in y dimension. denotes the magnetic flux through the tube, which leads to a phase factor in the inter-chain hopping integral as shown in Eq. (25). The relation of θ and is with being the magnetic flux quantum.

Taking the periodic boundary condition along y direction into account, the fermions are expressed in momentum space by the Fourier transformation

where k and q represent the momenta along x and y directions, respectively. Then taking thermodynamic limit while keeping M fixed, the Hamiltonian (25) is rewritten as
Here and the above summation is taken over q in the following set
The BdG Hamiltonian is
where
The form of the Hamiltonian in momentum space as shown in Eq. (27) obviously indicates that the system can be viewed as stack of M independent Kitaev chains with renormalized chemical potential and gap function as shown in Eqs. (30) and (31). Therefore by calculating the topological invariant of the BdG Hamiltonian for each given q, one can obtain the topological invariant of the whole tube system,

3.1. Class BDI

For and , the system as shown in Eq. (25) belongs to class BDI, which has time-reversal, particle-hole, and chiral symmetries.[34] Since , has chiral symmetry with for each given q. And thus can be transformed into off-diagonal form,

where
and . For each q, a winding number is associated with by
where . equals 1 if the following topological condition
is satisfied and 0 otherwise. Equation (36) shows that topological conditions for q and are the same and therefore . Correspondingly equation (32) becomes
For odd M, π is not in the set of q’s as shown in Eq. (28) and therefore there is no contribution to the topological invariant from as shown in Eq. (37).

The phase diagrams of the Kitaev tube models belonging to class BDI are illustrated in Fig. 4. The system is characterized by the Z index overall and the even–odd effect with respect to M is also manifested in the region where

similar to the Kitaev ladder belonging to the same class. Specifically speaking, for odd M the topological number is , while for even M the topological number is as dictated by Eq. (37).

Fig. 4. (color online) Phase diagrams of the Kitaev tubes composed of even and odd number of chains, which belongs to class BDI. The tube is composed of 6 Kitaev chains for the upper panel, while 5 chains for the lower panel. tx is chosen as energy unit. and . The integer numbers in the figure represent the index. In the regions where , the topological indices are 2,4,6 for M = 6 (upper panel) and 1,3,5 for M = 5 (lower panel).

Next we study the MZMs of the Kitaev tube system belonging to class BDI. For this purpose the Hamiltonian (25) is rewritten in the Majorana basis (15),

The MZM has the same form as Eq. (17) and the wave function also satisfies Eq. (18). However, we emphasize that here for the tube structure satisfies the periodic boundary condition in the y direction
instead of the open boundary condition as shown in Eq. (19) of the ladder structure. Therefore Fourier transformation is applied
and we have
with
Similar to the ladder case, we find the correspondence of bulk topological invariant and boundary zero modes in the Kitaev tube of class BDI. Iff certain q satisfies the topological condition Eq. (36) the solution of Eq. (41) exists (See Appendix A for detail) which ensures the presence of the boundary zero modes. If q equals 0 (or if M is even), can give rise to an MZM considering that is real. Otherwise, is complex because takes complex value. However since q and are counterparts of each other and , we have . From and , two complex fermions can be constructed according to and with relation . Consequently, two Majorana fermions and are obtained.

3.2. Class D

In this section, we consider the case or which results in that the system Hamiltonian Eq. (25) belongs to class D because the time-reversal and chiral symmetries are broken. For each given q which is not equal to the time-reversal-invariant momentum 0 or π, the topological invariant of the 1D BdG Hamiltonian is zero because it belongs to symmetry class A[34] as known from Eqs. (29), (30), and (31). For q = 0, π, the BdG Hamiltonian simply is

which belongs to class BDI and accordingly . Note that the Kitaev tube composed of odd number of chains has only one time-reversal-invariant momentum, i.e., q = 0, and thus the topological number is . Consequently there are two topologically distinct phases and the topological condition for is
On the other hand, the Kitaev tube composed of even number of chains has two time-reversal-invariant momenta, i.e., 0 and π, and thus the topological number is . Three topologically distinct phases exist and in addition to Eq. (44) for there is also a topological condition for
From the above results, the topological invariant of the class D Kitaev tube can be written as
As a result, the topological index is for odd M while 0,1,2 for even M.

The phase diagrams of the Kitaev tube models belonging to class D in the plane are plotted in Fig. 5. One can see that the phase diagram of even M has a period of π instead of the trivial period 2π for odd M. This finding can be explained by derived from Eq. (43). Correspondingly we have and thus for even M, which is compatible with the fact that the topological conditions (44) and (45) are interchanged through . From the relation and with B the magnetic field and S the sectional area of the tube, an oscillation of the topological properties as a function of the magnetic field B can be inferred with the period , which depends on the number of chains M and the sectional area S of the tube.

Fig. 5. (color online) Phase diagrams of the Kitaev tubes composed of even (a) and odd number (b) of chains, which belongs to class D. tx is chosen as energy unit. , , and . The integer numbers in panels (a) represent 0,1,2 index, while those in panel (b) represent the index.

At last we consider the MZMs of the Kitaev tube of class D. Two cases with odd and even chains are studied. In terms of Majorana basis, equation (38) reads

The results in momentum space suggest that the one solution of the MZM should be translation invariant along the y direction, namely corresponding to the time-reversal-invariant momentum q = 0 in Eq. (43). Therefore
and from , we have
The above equation is equivalent to that of the 1D Kitaev chain with an effective chemical potential
and there exists one MZM if the topological condition (44) is satisfied. (See Appendix A for detail)

If M is even, there may exist another MZM whose wave function is staggered along the y direction corresponding to the time-reversal-invariant momentum in Eq. (43), and accordingly

Substituting Eq. (51) into , we have Eq. (49) with
Accordingly the presence of MZM can be guaranteed iff the topological condition (45) is satisfied.

The bulk-boundary correspondence is reconfirmed in the Kitaev tube belonging to class D as in the cases of the Kitaev ladder and tube belonging to class BDI.

4. Conclusion

In this paper, we have studied two quasi-1D Kitaev models, namely the Kitaev ladder and tube, respectively. The topological phase diagrams and wave functions of the Majorana zero modes have been given explicitly. Under the condition that the amplitude of the interchain pairing potential is greater than that of the interchain hopping integral, we have revealed that the M-leg Kitaev ladder model of class BDI is characterized by the index if M is odd, while by the index if M is even. For the Kitaev tube model, the topological invariants have been given after taking Fourier transformation in the circumferential direction. The topological number of the Kitaev tube which belongs to class BDI also shows even–odd effect with respect to the number of chains in certain regions of the phase diagrams. For Kitaev tube of class D, the topological number is when the number of chains is odd, while 0, 1, 2 when the number of chains is even. For all the cases studied in this paper, we have shown that bulk topological invariant equals to the number of MZMs, which demonstrates the bulk-boundary correspondence in these systems.

Appendix A

In the study of MZMs, the Kitaev ladder or tube has a semi-infinite structure with an open boundary at the left-most lattice sites. As shown in the text of the paper, we find that the transformed wave function Fn satisfying the following second-order difference equation

Here , , and μ is generally a complex number. The general solution of Eq. (A1) is
where C is a constant to be determined by the normalization condition of the Majorana fermion and and are the roots of the following quadratic equation
with
The form of solution as Eq. (A2) obviously gurrentees the open boundary condition . Furthermore, as a boundary mode the MZM is localized around the left end of the semi-infinite ladder or tube, which demands that . To meet this, the following damping condition
should be fulfilled. After straightforward derivation one can easily find that the above damping condition is equivalent to
where and are the real and imaginary parts of μ, i.e. . As we show in the main text of the paper, the topological conditions obtained by examining the Hamiltonian in momentum space, i.e. Eqs. (12), (13), (36), (44), and (45), are the same as those deduced from the condition Eq. (A6) for the existence of MZMs in real space with the renormalized chemical potential μ being equal to Eqs. (23), (42), (50), and (52). This observation is exactly the manifestation of the bulk-boundary correspondence.

Reference
[1] Kitaev A Y 2001 Phys.-Usp. 44 131 http://stacks.iop.org/1063-7869/44/i=10S/a=S29
[2] Majorana E 1937 Nuovo Cimento 14 171 http://physics.princeton.edu/mcdonald/examples/EP/majorana_nc_14_171_37.pdf
[3] Kitaev A Y 2009 AIP Conf. Proc. 1134 22 http://aip.scitation.org/doi/abs/10.1063/1.3149495
[4] Alicea J 2012 Rep. Prog. Phys. 75 076501 http://stacks.iop.org/0034-4885/75/i=7/a=076501
[5] Fendley P 2012 J. Stat. Mech. 2012 11020 http://stacks.iop.org/1742-5468/2012/i=11/a=P11020
[6] Niu Y Chung S B Hsu C H Mandal I Raghu S Chakravarty S 2012 Phys. Rev. B 85 035110 https://link.aps.org/doi/10.1103/PhysRevB.85.035110
[7] DeGottardi W Thakurathi M Vishveshwara S Sen D 2013 Phys. Rev. B 88 165111 https://link.aps.org/doi/10.1103/PhysRevB.88.165111
[8] Greiter M Schnells V Thomale R 2014 Ann. Phys. 351 1026 http://www.sciencedirect.com/science/article/pii/S0003491614002498
[9] Hegde S S Shivamoggi V Vishveshwara S Sen D 2015 New J. Phys. 17 053036 http://stacks.iop.org/1367-2630/17/i=5/a=053036
[10] Hegde S S Vishveshwara S 2016 Phys. Rev. B 94 115166 https://link.aps.org/doi/10.1103/PhysRevB.94.115166
[11] Sau J D Clarke D J Tewari S 2011 Phys. Rev. B 84 094505 https://link.aps.org/doi/10.1103/PhysRevB.84.094505
[12] Alicea J Oreg Y Refael G von Oppen F Fisher M P A 2011 Nat. Phys. 7 412
[13] Leijnse M Flensberg K 2012 Semicond. Sci. Technol. 27 124003 http://stacks.iop.org/0268-1242/27/i=12/a=124003
[14] van Heck B Akhmerov A R Hassler F Burrello M Beenakker C W J 2012 New J. Phys. 14 035019 http://stacks.iop.org/1367-2630/14/i=3/a=035019
[15] Hyart T van Heck B Fulga I C Burrello M Akhmerov A R Beenakker C W J 2013 Phys. Rev. B 88 035121 https://link.aps.org/doi/10.1103/PhysRevB.88.035121
[16] Burnell F J 2014 Phys. Rev. B 89 224510 https://link.aps.org/doi/10.1103/PhysRevB.89.224510
[17] Fu L Kane C L 2008 Phys. Rev. Lett. 100 096407 https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.100.096407
[18] Oreg Y Refael G von Oppen F 2010 Phys. Rev. Lett. 105 177002 https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.105.177002
[19] Lutchyn R M Sau J D Das Sarma S 2010 Phys. Rev. Lett. 105 077001 https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.105.077001
[20] Stanescu T D Lutchyn R M Das Sarma S 2011 Phys. Rev. B 84 144522 https://link.aps.org/doi/10.1103/PhysRevB.84.144522
[21] Mourik V Zuo K Frolov S M Plissard S R Bakkers E P A M Kouwenhoven L P 2012 Science 336 1003 http://science.sciencemag.org/content/336/6084/1003
[22] Das A Ronen Y Most Y Oreg Y Heiblum M Shtrikman H 2012 Nat. Phys. 8 887 https://www.nature.com/articles/nphys2479.pdf
[23] Deng M T Yu C L Huang G Y Larsson M Caroff P Xu H Q 2012 Nano Lett. 12 6414 http://pubs.acs.org/doi/ipdf/10.1021/nl303758w
[24] Rokhinson L P Liu X Furdyna J K 2012 Nat. Phys. 8 795 https://www.nature.com/articles/nphys2429
[25] Potter A C Lee P A 2010 Phys. Rev. Lett. 105 227003 https://link.aps.org/doi/10.1103/PhysRevLett.105.227003
[26] Wimmer M Akhmerov A R Medvedyeva M V Tworzydlo J Beenakker C W J 2010 Phys. Rev. Lett. 105 046803 https://journals.aps.org/prl/abstract/10.1103/PhysRevLett.105.046803
[27] Zhou B Shen S Q 2011 Phys. Rev. B 84 054532 https://journals.aps.org/prb/pdf/10.1103/PhysRevB.84.054532
[28] DeGottardi W Sen D Vishveshwara S 2011 New J. Phys. 13 065028 http://stacks.iop.org/1367-2630/13/i=6/a=065028
[29] Wu N 2012 Phys. Lett. A 376 3530 https://ac.els-cdn.com/S0375960112010444/1-s2.0-S0375960112010444-main.pdf?_tid=40281c16-e861-11e7-9e4c-00000aab0f01&acdnat=1514089184_4afed96851cc84522cd9acd731c85548
[30] Wakatsuki R Ezawa M Nagaosa N 2014 Phys. Rev. B 89 174514 https://link.aps.org/doi/10.1103/PhysRevB.89.174514
[31] Wang H Shao L Pan Y Shen R Sheng L Xing D 2016 Phys. Lett. A 380 3936 http://www.sciencedirect.com/science/article/pii/S0375960116311628
[32] Nadj-Perge S Drozdov I K Li J Chen H Jeon S Seo J MacDonald A H Bernevig B A Yazdani A 2014 Science 346 602 http://science.sciencemag.org/content/346/6209/602.abstract
[33] Stanescu T D Tewari S 2013 J. Phys.: Condens. Matter 25 233201 http://iopscience.iop.org/article/10.1088/0953-8984/25/23/233201/pdf
[34] Chiu C K Teo J C Y Schnyder A P Ryu S 2016 Rev. Mod. Phys. 88 035005 https://link.aps.org/doi/10.1103/RevModPhys.88.035005
[35] Zhou B Z Zhou B 2016 Chin. Phys. B 25 107401 https://cpb.iphy.ac.cn/CN/abstract/article_68247.shtml